Reliability Estimate by Various Models

Author

Tzu-Yao Lin

Published

December 16, 2022

Code
#|label: setup
#|echo: false

knitr::opts_chunk$set(
  message = FALSE, warning = FALSE, collapse = TRUE
)

1 Dataset

This dataset is come from Koval and his colleagues’ (2013) study in which they use a novel paradigm and anlytical approach to study the dynamic process of depression. In past study, depression is associated with higher average levels of negative affect (NA) and lower average levels of positive affect (PA). However, few studies untangle the relationship between the pattern of affective fluctuations and the dynamics of depression. Therefore, in their study, they collected 99 (the final number) subjects’ depression symptom by the Center for Epidemiologic Studies Depression Scaleand (CES-D) and daily experiences of affect ratings 10 times/day for 7 days using the experience sampling method (ESM). Affect ratings were made on two items measuring PA (= mean(happy, relaxed)) and four items measuring NA (=mean(sad, depressed, anxious, angry)), with the scale from 1 (not at all) to 100 (very much).

Code
#|label: load-packages-and-rawdata

library(tidyverse)
library(lubridate)

rawdata <- read_tsv("data_1beep_no1st beep_annette.tsv")
rmarkdown::paged_table(rawdata)
Code
#|label: sorting-data

data <- rawdata %>% 
  mutate(Pos = PA, 
         Neg = `NA`,
         Date_Time = ymd_hms(str_glue("{year}-{month}-{day} {hour}:{min}:{sec}"), tz = "CET"),
         Date = as_date(Date_Time),
         Time = hms::as_hms(Date_Time),
         WDay = wday(Date, label = TRUE),
         Subject = factor(cumsum(PpID != lag(PpID, default = 0))), 
         .keep = "none") %>% 
  group_by(Subject) %>% 
  mutate(Day = factor(cumsum(Date != lag(Date, default = origin)))) %>% 
  group_by(Subject, Date) %>% 
  mutate(Moment = factor(1:n())) %>% 
  ungroup() %>% 
  pivot_longer(cols = c(Pos, Neg), names_to = "Affect", values_to = "Score") 

rmarkdown::paged_table(data)
Code
#|label: line-plot-of-score-in-diff-day-time 

data %>% 
  filter(Subject %in% c(6, 17, 44, 60)) %>%  
  drop_na() %>% 
  ggplot(aes(x = as_datetime(paste(today(), Time)),
             y = Score, 
             group = WDay, color = WDay)) + 
  geom_point() +
  geom_line() +
  scale_x_datetime(date_breaks = "2 hour", date_labels = "%H:%M", 
                   limits = c(as.POSIXct(c("9:30", "23:30"), format = "%H:%M"))) + 
  facet_grid(Subject ~ Affect)

Nevertheless, I found some discrepancies between the data I got and the descriptions in the literature.

  1. The number of subject in my data is 97;
  2. The sampling time points are not always 7 (days) and 10 (times). There are many missing values.

All analyses were conducted using the maximum available sample size: for analyses involving ESM data n=95; for all other analyses n=99.

Code
#|label: num-of-missing-data-1

t1 <- data %>% 
  group_by(Subject, Affect) %>% 
  summarise(num_day = length(unique(Day)),
            num_obser = n(),
            num_missing = sum(is.na(Score)),
            num_obser_valid = num_obser - num_missing) %>% 
  ungroup()
rmarkdown::paged_table(t1)
Code
summary(t1)
##     Subject       Affect             num_day        num_obser    
##  1      :  2   Length:194         Min.   :7.000   Min.   :62.00  
##  2      :  2   Class :character   1st Qu.:7.000   1st Qu.:66.00  
##  3      :  2   Mode  :character   Median :7.000   Median :66.00  
##  4      :  2                      Mean   :7.031   Mean   :66.75  
##  5      :  2                      3rd Qu.:7.000   3rd Qu.:67.00  
##  6      :  2                      Max.   :8.000   Max.   :78.00  
##  (Other):182                                                     
##   num_missing     num_obser_valid
##  Min.   : 0.000   Min.   :39.00  
##  1st Qu.: 3.000   1st Qu.:58.00  
##  Median : 6.000   Median :61.00  
##  Mean   : 6.737   Mean   :60.02  
##  3rd Qu.: 9.000   3rd Qu.:64.00  
##  Max.   :29.000   Max.   :73.00  
## 
table(t1$num_day)
## 
##   7   8 
## 188   6
Code
#|label: num-of-missing-data-2

t2 <- data %>% 
  group_by(Subject, Affect, Day) %>% 
  summarise(num_moment = n(), 
            num_missing = sum(is.na(Score)), 
            num_valid = num_moment - num_missing) %>% 
  ungroup()
rmarkdown::paged_table(t2)
Code
summary(t2)
##     Subject        Affect               Day        num_moment    
##  52     :  16   Length:1364        1      :194   Min.   : 5.000  
##  70     :  16   Class :character   2      :194   1st Qu.: 9.000  
##  87     :  16   Mode  :character   3      :194   Median :10.000  
##  1      :  14                      4      :194   Mean   : 9.494  
##  2      :  14                      5      :194   3rd Qu.:10.000  
##  3      :  14                      6      :194   Max.   :20.000  
##  (Other):1274                      (Other):200                   
##   num_missing        num_valid     
##  Min.   : 0.0000   Min.   : 0.000  
##  1st Qu.: 0.0000   1st Qu.: 8.000  
##  Median : 1.0000   Median : 9.000  
##  Mean   : 0.9582   Mean   : 8.536  
##  3rd Qu.: 1.0000   3rd Qu.:10.000  
##  Max.   :10.0000   Max.   :11.000  
## 
table(t2$num_moment)
## 
##   5   6   7   8   9  10  11  14  20 
##  10  20 138  42 188 864  98   2   2
table(t2$num_valid)
## 
##   0   2   3   4   5   6   7   8   9  10  11 
##   2   4   2  30  46  89 147 208 348 448  40

2 Generalizability Theory (GT)

2.1 Model

Moment is nested into Day

\[ N_{ik(j)} = \mu +S_i + D_j + M_{k(j)} + (SD)_{ij} + \epsilon_{ik(j)} \]

where

  • \(N_{ijk}\): the score of negative affect for subject \(i\) at moment \(k\) in day \(j\).

  • \(\mu\) : the grand mean.

  • \(S_i\): the random effect for the subject \(i\).

  • \(D_j\): the random effect for the day \(j\).

  • \(M_{k(j)}\): the random effect of the moment \(k\) and it is nested in the specific day \(j\).

The model in Schönbrodt et al. (2022) was designed as that Day and Moment are crossed each other. However, in our data, Moment nested within Day is more reasonable because the time points of each measurements were random in each day and there are no specific meaning of those sampling points.

2.2 Reliability

Follow the Schönbrodt et al. (2022) and Shrout & Lane (2012), the reliability of our model should be distinguished as between subject reliability and within subject reliability. However, these two reliabilities are defined by correlation format not by GT theory, which Schönbrodt et al. (2022) and Shrout & Lane (2012) approached (To be honest, I still cannot figure out how they built out their definition). Furthermore, in Schoenbrodt et al’s (2022) model, they assume that Moment is crossed with Day, the model has a Couple effect over subjects and item effect within each subject, and Moment effect and Item effect are fixed. So, You’ll notice a clear difference between the way I define these reliaiblities and the way they define them.

  1. Between subject reliability

\[ \bar{N}_{i} = \frac{1}{n_D n_{M_{ij}}}\sum_{j=1}^{n_D}\sum_{k=1}^{n_{M_{ij}}} N_{ik(j)}= \mu +S_i + \frac{1}{n_D}\sum_{j=1}^{n_D}D_j + \frac{1}{n_{M_{ij}}}\sum_{k=1}^{n_ {M_{ij}}} M_{k(j)} + \frac{1}{n_D}\sum_{j=1}^{n_D}(SD)_{ij} + \frac{1}{n_D n_{M_{ij}}}\sum_{j=1}^{n_D}\sum_{k=1}^{n_{M_{ij}}} \epsilon_{ik(j)} \]

\[\begin{align} R^{bw}_{abs} &= Corr(\bar{N}_{i}, \bar{N'}_{i}) \\ &= \frac{Cov \left( \begin{matrix} \mu +S_i + \frac{1}{n_D}\sum_{j=1}^{n_D}D_j + \frac{1}{n_{M_{ij}}}\sum_{k=1}^{n_ {M_{ij}}} M_{k(j)} + \frac{1}{n_D}\sum_{j=1}^{n_D}(SD)_{ij} + \frac{1}{n_D n_{M_{ij}}}\sum_{j=1}^{n_D}\sum_{k=1}^{n_{M_{ij}}} \epsilon_{ik(j)}, \\ \mu +S_i + \frac{1}{n_D}\sum_{j'=1}^{n_D}D_{j'} + \frac{1}{n_{M_{ij'}}}\sum_{k=1}^{n_ {M_{ij'}}} M_{k(j')} + \frac{1}{n_D}\sum_{j=1}^{n_D}(SD)_{ij'} + \frac{1}{n_D n_{M_{ij'}}}\sum_{j=1}^{n_D}\sum_{k=1}^{n_{M_{ij'}}} \epsilon_{ik(j')} \end{matrix} \right) }{\sqrt{Var(\bar{N}_{i})Var(\bar{N'}_{i})}} \\ &= \frac{\sigma^2_S}{\sigma^2_S + \sigma^2_D / n_D + \color{red}{\sigma^2_M} / n_M+ \sigma^2_{SD} / n_D + \sigma^2_{\epsilon} / n_D n_M} \end{align}\]

Warning

However, the numbers of moments for each subject are not same because of the missing data. Currently, I set all \(n_{M_{ij}}\)’s equal to 10 for convenience.

  1. Within subject change reliability

From day to day or from moment to moment

\[ \bar{N}_{ij} = \frac{1}{n_{M_{ij}}}\sum_{k=1}^{n_{M_{ij}}} N_{ik(j)}= \mu +S_i + D_j + \frac{1}{n_{M_{ij}}}\sum_{k=1}^{n_{M_{ij}}} M_{k(j)} + (SD)_{ij} + \frac{1}{n_{M_{ij}}}\sum_{k=1}^{n_{M_{ij}}} \epsilon_{ik(j)} \]

\[\begin{align} R^{WS}_{interday, abs} &= Corr(\bar{N}_{ij}, \bar{N}_{ij'} | S_i = s_i) \\ &= \frac{Cov \left( \begin{matrix} \mu +s_i + D_j + \frac{1}{n_{M_{ij}}}\sum_{k=1}^{n_{M_{ij}}} M_{k(j)} + (sD)_{ij} + \frac{1}{n_{M_{ij}}}\sum_{k=1}^{n_{M_{ij}}} \epsilon_{ik(j)}, \\ \mu +s_i + D_{j'} + \frac{1}{n_{M_{ij'}}}\sum_{k=1}^{n_{M_{ij'}}} M_{k(j')} + (sD)_{ij'} + \frac{1}{n_{M_{ij'}}}\sum_{k=1}^{n_{M_{ij}}} \epsilon_{ik(j')} \end{matrix} \right) }{\sqrt{Var(\bar{N}_{ij})Var(\bar{N}_{ij'})}}\\ &= \frac{\color{red}{??????}}{\sigma^2_D + \sigma^2_M/n_M + \sigma^2_{SD} + \sigma^2_{\epsilon}/n_M} \end{align}\]

\[\begin{align} R^{WS}_{intermoment, abs} &= Corr(Y_{ik(j)}, Y_{ik'(j)} | S_i = s_i) \\ &= \frac{Cov \left( \begin{matrix} \mu +s_i + D_j + M_{k(j)} + (sD)_{ij} + \epsilon_{ik(j)}, \\ \mu +s_i + D_{j} + M_{k'(j)} + (sD)_{ij} + \epsilon_{ik'(j)} \end{matrix} \right) }{\sqrt{Var(Y_{ik(j)})Var(Y_{ik'(j)})}}\\ &= \frac{\sigma^2_D + \sigma^2_{SD}}{\sigma^2_D + \sigma^2_M + \sigma^2_{SD} + \sigma^2_{\epsilon}} \end{align}\]

Compare to Schönbrodt et al. (2022)

Code
#|label: fitted-by-lme4

library(lme4)

neg_data <- data %>% 
  filter(Affect == "Neg",
         Day %in% 1:7) %>% 
  drop_na()

neg_lme <- lmer(Score ~ (1 | Subject) + (1 | WDay) + (1 | WDay:Moment) + (1 | Subject:WDay), 
                data = neg_data)
summary(neg_lme)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Score ~ (1 | Subject) + (1 | WDay) + (1 | WDay:Moment) + (1 |  
##     Subject:WDay)
##    Data: neg_data
## 
## REML criterion at convergence: 43793.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.1503 -0.4552 -0.1431  0.2612  5.8739 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  Subject:WDay (Intercept)  34.24    5.851  
##  Subject      (Intercept) 108.01   10.393  
##  WDay:Moment  (Intercept)   0.00    0.000  
##  WDay         (Intercept)   0.00    0.000  
##  Residual                  90.14    9.494  
## Number of obs: 5797, groups:  
## Subject:WDay, 678; Subject, 97; WDay:Moment, 80; WDay, 7
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   15.482      1.086   14.25
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

pos_data <- data %>% 
  filter(Affect == "Pos",
         Day %in% 1:7) 

pos_lme <- lmer(Score ~ (1 | Subject) + (1 | WDay) + (1 | WDay:Moment) + (1 | Subject:WDay), 
                data = pos_data)
summary(pos_lme)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Score ~ (1 | Subject) + (1 | WDay) + (1 | WDay:Moment) + (1 |  
##     Subject:WDay)
##    Data: pos_data
## 
## REML criterion at convergence: 49730.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0791 -0.5708  0.0766  0.6421  3.2019 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  Subject:WDay (Intercept)  57.993   7.615  
##  Subject      (Intercept) 155.594  12.474  
##  WDay:Moment  (Intercept)   3.788   1.946  
##  WDay         (Intercept)   4.509   2.123  
##  Residual                 262.216  16.193  
## Number of obs: 5794, groups:  
## Subject:WDay, 678; Subject, 97; WDay:Moment, 80; WDay, 7
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)    57.66       1.56   36.95
Warning

Caution: When I ran both model estimations, lme4 displayed a warming about “boundary (singular) fit”.

2.3 Result

Negative Affect Positive Affect
Between subject reliability \[ \frac{108.01}{108.01 + \frac{0}{7} + \frac{0}{10} + \frac{34.24}{7} + \frac{90.14}{7\times10}} = 0.946 \] \[ \frac{155.59}{155.59 + \frac{4.51}{7} + \frac{3.79}{10} + \frac{57.99}{7} + \frac{262.22}{7\times10}} = 0.923 \]
Within subject change reliability \[ \frac{0 + 34.24}{0 + 0 + 34.24 + 90.14} =0.275 \] \[ \frac{4.51 + 57.99}{4.51 + 3.79 + 57.99 + 262.22} = 0.190 \]
Code
# Not run yet
library(brms)

neg_brm <- brm(Score ~ (1 | Subject) + (1 | WDay) + (1 | WDay:Moment) + (1 | Subject:WDay),
               data = neg_data)
summary(neg_brm)

pos_brm <- brm(Score ~ (1 | Subject) + (1 | WDay) + (1 | WDay:Moment) + (1 | Subject:WDay),
               data = pos_data)
summary(pos_brm)

2.4 Changing variables

Code
ggplot(data, aes(x = as_datetime(paste(today(), Time)))) + 
  geom_histogram(color = "white", binwidth = 30*60, boundary = 0) +
  scale_x_datetime(name = "Time", date_breaks = "2 hours", date_labels = "%H:%M") 

  • morning: 9:00-12:00

  • early-afternoon: 12:00-15:00

  • late-afternoon: 15:00-18:00

  • early-evening: 18:00-21:00

  • late-evening: 21:00-24:00

Code
data2 <- data %>% 
  filter(Day %in% 1:7) %>% 
  mutate(Weekend = WDay %in% c("Sat", "Sun"), 
         CMoment = factor(cut(as.double(Time), 
                              breaks = c(9, 12, 15, 18, 21, 24) * 3600, 
                              label = FALSE),
                          labels = c("morning", 
                                     "early-afternoon", "late-afternoon",
                                     "early-evening", "late-evening"),
                          order = TRUE)) %>% 
  group_by(Subject, Weekend, Affect, CMoment) %>% 
  summarise(N = sum(!is.na(Score)),
            AveScore = mean(Score, na.rm = TRUE))
   
rmarkdown::paged_table(data2)
Code

data2 %>% filter(Affect =="Neg") %>% 
  group_by(CMoment) %>% 
  summarise(n_ = n(), 
            N = sum(N))
## # A tibble: 5 × 3
##   CMoment            n_     N
##   <ord>           <int> <int>
## 1 morning           189   853
## 2 early-afternoon   194  1375
## 3 late-afternoon    194  1515
## 4 early-evening     194  1455
## 5 late-evening      189   599

2.4.1 Model

\[ N^*_{ik(j)} = \mu +S_i + W_j + M^*_{k(j)} + (SW)_{ij} + \epsilon_{ik(j)} \]

where

  • \(N^*_{ijk}\): the new average score of negative affect for subject \(i\) at moment \(k\) in weekend \(j=1\) or not \(j=0\)

  • \(\mu\) : the grand mean.

  • \(S_i\): the random effect for the subject \(i\).

  • \(W_j\): the fixed effect to indicate whether the observation is in weekend or not.

  • \(M_{k(j)}\): the random effect of the new categorical moment $k$, including: morning, early-afternoon, late-afternoon, early-evening, and late-evening, with 5 levels. Moment* is nested in \(W_j\).

2.4.2 Reliability

  1. Between subject reliability

\[ \begin{align}R^{bw}_{abs} &= Corr(\bar{N}^*_{i}, \bar{N'}^*_{i} \mid W_j = w_j) \\&= \frac{Cov \left( \begin{matrix} \mu +S_i + \frac{1}{n_W}\sum_{j=1}^{n_W}w_j + \frac{1}{n_{M^*_{ij}}}\sum_{k=1}^{n_{M^*_{ij}}} M^*_{k(j)} + \frac{1}{n_W}\sum_{j=1}^{n_W}(Sw)_{ij} + \frac{1}{n_W n_{M^*_{ij}}}\sum_{j=1}^{n_D}\sum_{k=1}^{n_{M^*_{ij}}} \epsilon_{ik(j)}, \\ \mu +S_i + \frac{1}{n_W}\sum_{j=1}^{n_W}w_{j} + \frac{1}{n_{M^*_{ij}}}\sum_{k=1}^{n_{M^*_{ij}}} M^*_{k'(j)} + \frac{1}{n_W}\sum_{j=1}^{n_W}(Sw)_{ij} + \frac{1}{n_W n_{M^*_{ij'}}}\sum_{j=1}^{n_W}\sum_{k=1}^{n_{M^*_{ij'}}} \epsilon_{ik'(j)} \end{matrix} \right) }{\sqrt{Var(\bar{N}^*_{i})Var(\bar{N'}^*_{i})}} \\&= \frac{\sigma^2_S + \sigma^2_{SW}/n_{SW}}{\sigma^2_S + \color{red}{\sigma^2_M} / n_M+ \sigma^2_{SW} / n_W + \sigma^2_{\epsilon} / n_W n_M}\end{align} \]

  1. Within subject change reliability

Because Weekend is a fixed effect, I’m not sure how to build up reliability under this situation.

2.4.3 Result

Code
neg_data2 <- data2 %>% filter(Affect == "Neg")
neg_lme2 <- lmer(AveScore ~ 1 + Weekend + (1 | Subject) + (1 | Weekend:CMoment) + (1 | Subject:Weekend), 
                data = neg_data2)
summary(neg_lme2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: AveScore ~ 1 + Weekend + (1 | Subject) + (1 | Weekend:CMoment) +  
##     (1 | Subject:Weekend)
##    Data: neg_data2
## 
## REML criterion at convergence: 6415
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1753 -0.4418 -0.1000  0.3452  5.7003 
## 
## Random effects:
##  Groups          Name        Variance Std.Dev.
##  Subject:Weekend (Intercept)  14.99    3.872  
##  Subject         (Intercept) 104.09   10.203  
##  Weekend:CMoment (Intercept)   0.00    0.000  
##  Residual                     31.88    5.646  
## Number of obs: 946, groups:  
## Subject:Weekend, 194; Subject, 97; Weekend:CMoment, 10
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  15.3918     1.1375  13.531
## WeekendTRUE  -0.1508     0.6671  -0.226
## 
## Correlation of Fixed Effects:
##             (Intr)
## WeekendTRUE -0.291
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

pos_data2 <- data2 %>% filter(Affect == "Pos")
pos_lme2 <- lmer(AveScore ~ 1 + Weekend + (1 | Subject) + (1 | Weekend:CMoment) + (1 | Subject:Weekend), 
                data = pos_data2)
summary(pos_lme2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: AveScore ~ 1 + Weekend + (1 | Subject) + (1 | Weekend:CMoment) +  
##     (1 | Subject:Weekend)
##    Data: pos_data2
## 
## REML criterion at convergence: 7355.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0769 -0.5248  0.0424  0.5137  3.5369 
## 
## Random effects:
##  Groups          Name        Variance Std.Dev.
##  Subject:Weekend (Intercept)  17.325   4.162  
##  Subject         (Intercept) 148.001  12.166  
##  Weekend:CMoment (Intercept)   6.134   2.477  
##  Residual                     96.931   9.845  
## Number of obs: 946, groups:  
## Subject:Weekend, 194; Subject, 97; Weekend:CMoment, 10
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   57.175      1.770  32.305
## WeekendTRUE    1.669      1.795   0.929
## 
## Correlation of Fixed Effects:
##             (Intr)
## WeekendTRUE -0.506
Negative Affect Positive Affect
Between subject reliability \[ \frac{104.09 + \frac{14.99}{2}}{104.09 + \frac{14.99}{2} + \frac{0}{5} + \frac{31.88}{2\times5}} = 0.972 \] \[ \frac{148.00 + \frac{17.33}{2}}{148.00 + \frac{17.33}{2} + \frac{6.13}{5} + \frac{96.93}{2\times5}} = 0.935 \]

3 Multistate-Singletrait Model (MSST)

3.1 Model

\[\begin{align} N_{it} &= \tau_{it} + \epsilon_{it} \\ &= \mu_{it} + \lambda_{T_{it}} \xi + \lambda_{S_{it}} \zeta_t + \epsilon_{it} \end{align}\]

  • \(N_{it}\): the score of negative affect for subject \(i\) at date-time \(t\).
  • \(\mu_{it}\): the intercept (the ground mean).
  • \(\xi\): the factor score of the latent trait variable (the stable part).
  • \(\zeta_t\): the factor score of the latent state variable (the deviate part).
  • \(\lambda_{T_{it}}\), \(\lambda_{S_{it}}\): the factor loadings for the the the latent Trait variable and latent State variable, respectively.
  • \(\epsilon_{it} \overset{i.i.d.}{\sim} N(0, \sigma^2_\epsilon)\): the random measurement error.

3.2 Reliability

  1. Reliability: \(Rel = 1 - \frac{Var(\epsilon_{it})}{Var(N_{it})}\)
  2. Consistency: \(Con = \frac{\lambda^2_{T_{it}}Var(\xi)}{Var(N_{it})}\)
  3. Occasion-specificity: \(Spe = \frac{\lambda^2_{S_jt}Var(\zeta_{t})}{Var(N_{it})}\)

4 Common-Unique Trait–State Model (CUTS)

4.1 Model

\[\begin{align} N_{it} &= T_{i} + S_{it} \\ &= (\mu_{it} + \lambda_{T_{it}} \xi + \vartheta_i) + (\lambda_{S_{it}} \zeta_t + \epsilon_{it}) \end{align}\]

  • \(N_{it}\): the score of negative affect for subject \(i\) at date-time \(t\).
  • \(T_i\): the trait variable
    • \(\mu_{i}\): the intercept (the ground mean).

    • \(\xi\): the common trait score.

    • \(\vartheta_j\): the unique trait score for \(i\).

  • \(S_{it}\): the state variable
    • \(\zeta_t\): the common state score.

    • \(\epsilon_{it} \overset{i.i.d.}{\sim} N(0, \sigma^2_\epsilon)\): the unique state score.

4.2 Reliability

  1. Reliability: \(Rel = 1 - \frac{Var(\epsilon_{it})}{Var(N_{it})}\)
  2. Total consistency: \(TCon = \frac{\lambda^2_{T_{it}}Var(\xi) + \vartheta_i}{Var(N_{it})}\)
    1. Common consistency: \(CCon = \frac{\lambda^2_{T_{it}}Var(\xi)}{Var(N_{it})}\)

    2. Unique consistency: \(UCon = \frac{Var(\vartheta_i)}{Var(N_{it})}\)

  3. Occasion-specificity: \(Spe = \frac{\lambda^2_{S_{it}}Var(\zeta_t)}{Var(N_{it})}\)

5 Trait–State-Occasion Model (TSO)

5.1 Model

\[\begin{align} N_{it} &= \tau_{it} + \epsilon_{it} \\ &= \mu_{it} + \lambda_{T_{it}}\xi_i + \lambda_{S_{jt}}O_t + \epsilon_{it} \\ O_t &= \beta_{t-1, t}O_{t-1} + \zeta_t \\ O_1 &= \zeta_1 \end{align}\]

  • \(N_{it}\): the score of negative affect for subject \(i\) at date-time \(t\).
  • \(\mu_{it}\): the intercept (the ground mean).
  • \(\xi_i\): the factor score of the latent trait variable of subject \(i\) (the stable part).
  • \(O_t\): the latent occasion-specific variable at time \(t\) (the deviate part which can be predictable).
  • \(\zeta_t\): the latent occasion-specific residual at time \(t\) (the deviate part but not be predictable).
  • \(\beta_{t-1, t}\): the auto-regressive effect between two consecutive occasions, \(t-1\) and \(t\).
  • \(\lambda_{T_{it}}\), \(\lambda_{S_{it}}\): the factor loadings for the the the latent Trait variable and latent State variable, respectively.
  • \(\epsilon_{it} \overset{i.i.d.}{\sim} N(0, \sigma^2_\epsilon)\): the random measurement error.

5.2 Reliability

  1. Reliability:\(Rel = 1 - \frac{Var(\epsilon_{it})}{Var(N_{it})}\)
  2. Consistency: \(Con = \frac{\lambda^2_{T_{it}Var(\xi_t)} + \beta^2_{t-1, t}\lambda^2_{S_it}Var(O_{t-1})}{Var(N_{it})}\)
    1. Predictability by trait: \(Con = \frac{\lambda^2_{T_{it}Var(\xi_t)}}{Var(N_{it})}\)

    2. Unpredictability by trait: \(Con = \frac{\beta^2_{t-1, t}\lambda^2_{S_it}Var(O_{t-1})}{Var(N_{it})}\)

  3. Occasion-specificity: \(Spe = \frac{\lambda^2_{S_{it}}Var(\zeta_t)}{Var(N_{it})}\)
Note

I found the three kinds of variance coefficients are hard to use correlations to define them.

6 Measurement Error Vector Autoregressive with Order 1 Model (MEVAR(1))

6.1 Model

Level 1 (within subject)

  • Measurement equation

\[\begin{pmatrix}P_{it} \\ N_{it}\end{pmatrix} = \begin{pmatrix} y_{1it} \\ y_{2it}\end{pmatrix} = \begin{pmatrix}\mu_{1i} \\ \mu_{2i}\end{pmatrix} + \begin{pmatrix}\tilde{y}_{1it} \\ \tilde{y}_{2it}\end{pmatrix} + \begin{pmatrix}\epsilon_{1it} \\ \epsilon_{2it}\end{pmatrix} \]

\[ \begin{pmatrix}\epsilon_{1it} \\ \epsilon_{2it}\end{pmatrix} \sim N_2\left(\begin{pmatrix}0 \\ 0\end{pmatrix}, \begin{pmatrix} \sigma^2_{\epsilon 11i} & \sigma^2_{\epsilon 12i} \\ \sigma^2_{\epsilon 12i} & \sigma^2_{\epsilon 22i}\end{pmatrix}\right) \]

  • Transition equation (state space model representation)

\[ \begin{pmatrix}\tilde{y}_{1it} \\ \tilde{y}_{2it}\end{pmatrix} = \begin{pmatrix} \phi_{11i} & \phi_{12i} \\ \phi_{12i} & \phi_{22i} \end{pmatrix} \begin{pmatrix}\tilde{y}_{1it-1} \\ \tilde{y}_{2it-1}\end{pmatrix} + \begin{pmatrix}\omega_{1it} \\ \omega_{2it}\end{pmatrix} \]

\[ \begin{pmatrix}\omega_{1it} \\ \omega_{2it}\end{pmatrix} \sim N_2 \left(\begin{pmatrix}0 \\ 0\end{pmatrix}, \begin{pmatrix} \sigma^2_{\omega 11i} & \sigma^2_{\omega 12i} \\ \sigma^2_{\omega 12i} & \sigma^2_{\omega 22i}\end{pmatrix}\right) \]

Level 2 (between subject)

\[ \left(\begin{array}{c}\mu_{1 i} \\\mu_{2 i} \\\Phi_{11 \mathrm{i}} \\\Phi_{12 \mathrm{i}} \\\Phi_{21 \mathrm{i}} \\\Phi_{22 \mathrm{i}}\end{array}\right) \sim M v N\left(\left(\begin{array}{c}\gamma_{\mu 1} \\\gamma_{\mu 1} \\\gamma_{\Phi 11} \\\gamma_{\Phi 12} \\\gamma_{\Phi 21} \\\gamma_{\Phi 22}\end{array}\right), \boldsymbol{\Psi} =\begin{pmatrix} \underset{2\times2}{\boldsymbol{\Psi}_\mu} & \underset{2\times4}{\boldsymbol{\Psi}_{\mu\Phi}} \\ \underset{4\times2}{\boldsymbol{\Psi}^\top_{\mu\Phi}} & \underset{4\times4}{\boldsymbol{\Psi}_\phi}\end{pmatrix}\right) \]

\[ \left(\begin{array}{l}\sigma_{\omega_{11 i}}^{2} \\\sigma_{\omega_{12 i}}^{2} \\\sigma_{\omega_{22 i}}^{2}\end{array}\right) \sim M v N\left(\left(\begin{array}{l}\gamma_{\sigma_{\omega 11}^{2}} \\\gamma_{\sigma_{\omega 12}^{2}} \\\gamma_{\sigma_{\omega 22}^{2}}^{2}\end{array}\right), \boldsymbol{\Psi}_{\sigma_{\omega}^{2}}\right) \]

\[ \left(\begin{array}{l}\sigma_{\epsilon_{11 i}}^{2} \\\sigma_{\epsilon_{12 i}^{2 i}}^{2} \\\sigma_{\epsilon_{22 i}}^{2}\end{array}\right) \sim M v N\left(\left(\begin{array}{l}\gamma_{\sigma_{\epsilon 11}^{2}} \\\gamma_{\sigma_{\epsilon 12}^{2}}^{2} \\\gamma_{\sigma_{\epsilon 22}^{2}}^{2}\end{array}\right), \boldsymbol{\Psi}_{\boldsymbol{\sigma}_{\epsilon}^{2}}\right) \]

6.2 Reliability

Total variance of data: \(Var(\mathbf{y}_{ti}) = \boldsymbol{\Psi}^2_\mu + E_i[\mathbf{T}_{\tilde{y}i}] +\boldsymbol{\gamma}_{\sigma^2_\epsilon}\)

Proof: (\(\mathbf{y}_{ti} = \begin{pmatrix} y_{1ti} \\ y_{2ti} \end{pmatrix}\) vector representation)

\[ \begin{align*}Var(\mathbf{y}_{ti}) &= Var(\boldsymbol{\mu}_i + \mathbf{\tilde{y}}_{ti} + \boldsymbol{\epsilon}_{ti} ) \\&= Var(\boldsymbol{\mu}_i) + Var(\mathbf{\tilde{y}}_{ti}) + Var(\boldsymbol{\epsilon}_{ti}) ) \\&= Part 1 + Part 2 + Part 3\end{align*} \]

Part 1:

\[\begin{align*}Var(\boldsymbol{\mu}_i)&= Var_i(E_t[\boldsymbol{\mu}_i \mid I=i]) + E_i[Var_t(\boldsymbol{\mu}_i \mid I=i)] \\&= Var_i(\boldsymbol{\mu}_i) + E_i[0] \\&= \begin{pmatrix} \psi^2_{\mu11} \psi_{\mu12} \\ \psi_{\mu12} \psi^2_{\mu22} \end{pmatrix} = \boldsymbol{\Psi}^2_\mu \end{align*}\]

Part 2:

\[ \mathbf{\tilde{y}}_{ti} = \boldsymbol{\Phi}\mathbf{\tilde{y}}_{t-1 i} + \boldsymbol{\omega}_{ti} = \boldsymbol{\Phi}^{t}\mathbf{\tilde{y}}_{0 i} + \sum_{k=1}^{t}\boldsymbol{\Phi}^{t-1}\boldsymbol{\omega}_{ti} \\Var_t(\mathbf{\tilde{y}}_{ti}) = \boldsymbol{\Phi}Var_t(\mathbf{\tilde{y}}_{t-1 i})\boldsymbol{\Phi}^\top + Var_t(\boldsymbol{\omega}_{ti}) \]

Consider \(t\) is at steady state. Let \(\mathbf{T}_{\tilde{y}i} := Var_t(\mathbf{\tilde{y}}_{ti})\), (cf., Kim & Nelson,1999, p. 27)

\[\begin{align*} \mathbf{T}_{\tilde{y}i} &= \boldsymbol{\Phi}\mathbf{T}_{\tilde{y}i}\boldsymbol{\Phi}^\top + \boldsymbol{\Sigma_{\omega i}} \\vec(\mathbf{T}_{\tilde{y}i}) &= vec(\boldsymbol{\Phi}\mathbf{T}_{\tilde{y}i}\boldsymbol{\Phi}^\top) + vec(\boldsymbol{\Sigma_{\omega i}}) \\ vec(\mathbf{T}_{\tilde{y}i}) &= (\boldsymbol{\Phi} \otimes \boldsymbol{\Phi}) vec(\mathbf{T}_{\tilde{y}i})+ vec(\boldsymbol{\Sigma_{\omega i}}) \quad (\text{by Lemma 1} \\vec(\mathbf{T}_{\tilde{y}i}) &= (\mathbf{I}-(\boldsymbol{\Phi} \otimes \boldsymbol{\Phi})^{-1} )vec(\boldsymbol{\Sigma_{\omega i}}) \\\mathbf{T}_{\tilde{y}i} &= mat((\mathbf{I}-(\boldsymbol{\Phi} \otimes \boldsymbol{\Phi})^{-1} )vec(\boldsymbol{\Sigma_{\omega i}})) \end{align*}\]

Than, we have

\[\begin{align*} Var(\mathbf{\tilde{y}}_{ti})&= Var_i(E_t[\mathbf{\tilde{y}}_{ti} \mid I=i]) + E_i[Var_t(\mathbf{\tilde{y}}_{ti}\mid I=i)] \\&= Var_i(\boldsymbol{\Phi}^{t}\mathbf{\tilde{y}}_{0 i}) + E_i[\mathbf{T}_{\tilde{y}i}] \\ &= E_i[\mathbf{T}_{\tilde{y}i}] \end{align*}\]

Part 3:

\[\begin{align*} Var(\boldsymbol{\epsilon}_{ti}))&= Var_i(E_t[\boldsymbol{\epsilon}_{ti} \mid I=i]) + E_i[Var_t(\boldsymbol{\epsilon}_{ti} \mid I=i)] \\&= Var_i(\mathbf{0}) + E_i \left[\begin{pmatrix} \sigma^2_{\epsilon11i} \sigma^2_{\epsilon12i} \\ \sigma^2_{\epsilon12i} \sigma^2_{\epsilon22i}\end{pmatrix}\right] \\ &= \begin{pmatrix} \gamma_{\sigma^2_{\epsilon11}} \gamma_{\sigma^2_{\epsilon12}} \\ \gamma_{\sigma^2_{\epsilon12}} \gamma_{\sigma^2_{\epsilon22}} \end{pmatrix} = \boldsymbol{\gamma}_{\sigma^2_\epsilon} \end{align*}\]

  • Reliability for systematic between-subject difference

\[\begin{align*} R_B(\mathbf{y}) &= \frac{\Psi^2_\mu}{Var(\mathbf{y})} \\ &=\frac{Cov\begin{pmatrix} \boldsymbol{\mu_i} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{it-1} + \boldsymbol{\omega}_{it} + \boldsymbol{\epsilon}_{it}, \\ \boldsymbol{\mu_{i'}} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{i't-1} + \boldsymbol{\omega}_{i't} + \boldsymbol{\epsilon}_{i't} \end{pmatrix}}{\sqrt{Var(\boldsymbol{\mu_i} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{it-1} + \boldsymbol{\omega}_{it} + \boldsymbol{\epsilon}_{it})Var(\boldsymbol{\mu_{i'}} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{i't-1} + \boldsymbol{\omega}_{i't} + \boldsymbol{\epsilon}_{i't})}} \\ &\color{blue}{= Corr(\mathbf{y}_{it}, \mathbf{y}_{i't})} \end{align*}\]

  • Reliability for within-subject fluctuations

\[\begin{align*} R_w(y_i) &= \frac{\tau^2_{\tilde{y}i}}{v(y_i)} \\ &= \frac{???}{\boldsymbol{\Phi}Var_t(\mathbf{\tilde{y}}_{t-1 i})\boldsymbol{\Phi}^\top + Var_t(\boldsymbol{\omega}_{ti}) + \sigma^2_{\epsilon_i}} \\ &=\frac{Cov\left(\begin{matrix} \boldsymbol{\mu_i} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{it-1} + \boldsymbol{\omega}_{it} + \boldsymbol{\epsilon}_{it}, \\ \boldsymbol{\mu_i} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{it'-1} + \boldsymbol{\omega}_{it'} + \boldsymbol{\epsilon}_{it'} \end{matrix} \mid i \right)}{\sqrt{Var(\boldsymbol{\mu_i} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{it-1} + \boldsymbol{\omega}_{it} + \boldsymbol{\epsilon}_{it} i)Var(\boldsymbol{\mu_i} + \boldsymbol{\Phi}\mathbf{\bar{y}}_{it'-1} + \boldsymbol{\omega}_{it'} + \boldsymbol{\epsilon}_{it'} \mid i )}} \\ &\color{blue}{= Corr(\mathbf{y}_{it}, \mathbf{y}_{it'} \mid i)} \end{align*}\]

where \(v(y_i) = \tau_{\tilde{y}i} + \sigma^2_{\epsilon i}\) and \(\tau_{\tilde{y}i}\) is the person-specific variance for variable \(\tilde{y}^j\) and it is equal to the \(j\)th diagonal element for that variable of the covariance matrix \(\mathbf{T}_{\tilde{y}i}\) for person i.

References

Koval, P., Pe, M. L., Meers, K., & Kuppens, P. (2013). Affect dynamics in relation to depressive symptoms: Variable, unstable or inert? Emotion, 13(6), 1132–1141. https://doi.org/10.1037/a0033579
Schönbrodt, F. D., Zygar-Hoffmann, C., Nestler, S., Pusch, S., & Hagemeyer, B. (2022). Measuring motivational relationship processes in experience sampling: A reliability model for moments, days, and persons nested in couples. Behavior Research Methods, 54(4), 1869–1888. https://doi.org/10.3758/s13428-021-01701-7
Shrout, P. E., & Lane, S. P. (2012). Reliability. 643–660. https://doi.org/10.1037/13619-034